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In this paper, we present a real non-linear differential equation for the two flavor neutrino oscil- 
lation problem in matter with an arbitrary density profile. We also present an exact series solution 
to this non-linear differential equation. In addition, we investigate numerically the convergence 



■ of this solution for different matter density profiles such as constant and linear profiles as well as 
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the Preliminary Reference Earth Model describing the Earth's matter density profile. Finally, we 
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' discuss other methods used for solving the neutrino flavor evolution problem. 
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I. INTRODUCTION 



In general, there are several phenomena and processes in physics, but also in other fields 
of science such as chemistry, that can be described in terms of a system with two (quantum 
mechanical) states and a time-dependent Hamiltonian, i.e., so-called two-level systems - 
neutrino oscillations with two flavors being one such system. Other representatives of such 
systems are, for example, a spin 1/2 particle in a time-dependent electromagnetic field 
having the states 'spin up' and 'spin down', K°-K° mixing, a Josephson device, nuclear 
magnetic resonance used for encoding bits of information (i.e., quantum bits for a quantum 
computer), the left and right chirality states of molecules in chemistry, etc. The problem of 
neutrino oscillations in matter, which we are concerned with in this paper, is mathematically 
equivalent to a spin 1/2 particle in a magnetic field that is constant in one direction, zero 
in another direction, and time-dependent in the last direction 111. 

Neu tri „„ orations have — been e Xt en 8i ve. y 8tudled in the h— flflflflf | 
and they act as the most plausible description of both the solar Q] and atmospheric 2| 
neutrino problems. At an early stage, neutrino oscillations were mainly investigated with 
two flavors and without including matter effects. Nowadays, we know that there are at 
least three neutrino flavors and that matter effects are important. For example, in matter, 
the so-called Mikheyev-Smirnov-Wolfenstein (MSW) effect j?l Q can take place, which is 
an amplifying resonant effect due to the presence of matter. However, in most situations, 
neutrino oscillations can be effectively investigated with two flavors, since the leptonic mixing 
in the 1-3 sector is indeed small leading to the fact that the full three flavor scenario can 
be decoupled into two effective two flavor scenarios, each of which can be studied separately. 

In this paper, we present an exact analytic solution to the two flavor neutrino oscillation 
problem in matter. Since there are many similar two-level systems, as discussed above, our 
solution will also be interesting and applicable to this kind of systems. However, before 
we proceed to present our solution, we will give a brief overview of what has previously 
been done in this field. Note that this overview is not presented in chronological order. 
First, in Refs. [9|, 110|], the neutrino flavor evolution has been investigated by a discretization 
of the effective potential. Second, exact solutions exist for a number of specific effective 



potentials 



ctive potential, becona, exact somt 
Q, H 111 3- Third, in Refs - Q 



151 ] . the evolution was studied by using an 



adiabatic approximation. Fourth, approximate solutions valid for small effective potentials 
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have recently been studied in detail jig . Finally, there have also been other attempts to 
write the evolution in terms of a second order non-linear ordinary differential equation jl7|. 
However, this has been done for the neutrino oscillation probability amplitudes and not for 
the neutrino oscillation probabilities. The advantage of working with a non-linear differential 
equation for the oscillation probability rather than a linear system of differential equations 
for the probability amplitudes is that we only have one real variable instead of two complex 
variables. The disadvantage is that the resulting differential equation is non-linear. 

This paper is organized as follows. In Sec. |Hj the neutrino flavor evolution in matter with 
two flavors is studied and a second order non-linear ordinary differential equation for the 
neutrino oscillation probability is derived. Then, in Sec. Illl| we perform series expansions 
of both the neutrino oscillation probability and the effective potential in order to solve 
the differential equation presented in Sec. HTJ Next, in Sec. IIVI we continue by studying 
the numerical convergence of the solution for a number of different effective potentials and 
baselines. In Sec. EJ we present a brief summary of other methods for solving the neutrino 
evolution in matter. Finally, in Sec. IVI| we summarize our results and give our conclusions. 



II. NEUTRINO FLAVOR EVOLUTION IN MATTER 

When neutrinos propagate in matter, neutrino flavors are affected differently by coherent 
forward scattering against the matter constituents. Assuming that there are no sterile 
neutrinos, the effect of matter is to add an effective potential to z/ e , this effective potential 
is given by V(t) = ^/2GpN e (t), where Gp is the Fermi coupling constant and N e (t) is the 
electron number density. 

In the two flavor case, the time evolution of a neutrino state \v{t)) = (\v e (t)) |z/j.(t))) T is 
given by 

i^ffil = (H yae + H aBt (t))\v(t)), (1) 

where if vac = ?7diag(m^, m^jU^ /2E is the free Hamiltonian in vacuum, H ma , t (t) = 
diag(V(t), 0) is the addition to the free Hamiltonian due to matter effects, and 

-(::) 

is the leptonic mixing matrix in vacuum. Here c = cos 8, s = sin 6, and 9 is the leptonic 
mixing angle. Adding or subtracting terms proportional to the unity operator to the total 



Hamiltonian H(t) = H vabC + Hmstit) will only contribute with an overall phase to the neutrino 
state \v{t)), and thus, does not affect the neutrino oscillation probabilities. Using this fact, 
the total Hamiltonian may be written as 



H(t) 



sin 26 



1 / V(t) - ^£ cos 26 
2^ ^fsin2# ^fcos26-V(t) 
1 " 



2E 2E 
Am 2 ( Am^ 

(Ji^rr sin 26 + a 3 V(t) - -^cos2# 



■ie ' \ - - it : ^ 

where the a^s (i = 1, 2, 3) are the Pauli matrices and Am 2 = m\ — m\ is the mass squared 
difference between the two mass eigenstates in vacuum. 

The density matrix pit) = \u{t)) {u{t)\ can be parameterized as pit) = (1 + S{t) • cr)/2, 
where 1 is the unity matrix, cr = (ctl a 2 o" 3 ) T is the vector of Pauli matrices, and S(t) is a 
vector such that S(t) 2 = 1. Differentiating the density matrix p with respect to time t, the 
equation of motion for S(t) becomes 



S(t) = S(t) x B(t), 



(4) 



where B(t) = ge x + f(t)e 3 and we have defined g = — sin(2^)Am 2 /2_E and f(t) = 
cos(26)Am 2 /2E — V(t). Note that 5 is independent of time t. The probability of neu- 
trinos produced as v e to oscillate into u x (where u x is some linear combination of and v T ) 
is now given by P(u e — > i/ x ) = P ea , = (1 — S s (t))/2. With the parameterization 



^ sin ct cos (3 ^ 



sin a sin /3 
cos a 



(5) 



where a = a(t) and (3 = Pit), we obtain the following non-linear system of ordinary differ- 
ential equations 



$ — g cot a cos (3 — f, 
a = g sin (3. 

Eliminating (3 from the above expressions, we obtain the differential equation 

[a + cot a{a 2 - G)] 2 = F(t)(G - a 2 ), 



(6) 
(7) 

(8) 



where F(t) = fit) 2 and G = g 2 . 



Now, we make the substitution p = S^it) = cos a after which Eq. (JHJ) becomes 



(p + Gp) 2 = F{t)[G{l-p 2 )-p 2 }. 



(9) 



Note that F = corresponds to the so-called MSW resonance condition cos(2#)Am 2 jlE = 
V. In this case, Eq. © takes the simple form 



with the trivial solutions p = Acos(gt) + Bsin(gt) just as expected. 

In general, the expression for P ex is known for constant matter density and is given by 



where 9 is the effective leptonic mixing angle in matter and Am 2 is the effective mass squared 
difference in matter. Using that p = 1 — 2P ex , it is a matter of trivial computation to show 
that this is the solution to Eq. (jOJ) with constant F, which corresponds to any constant 
matter density. 

In the three (n) flavor case, the density matrix can be parameterized by four [2(n — 1)] 
real parameters. If we would adopt our approach to the three (n) flavor case, then we would 
end up with a system of seven [2(n — 1)] non-linear ordinary differential equations, which, 
in principle, can be solved in a manner analogous to the one described above for the two 
flavor case. 

III. SERIES EXPANSION OF THE SOLUTION 

In order to solve the propagation of neutrinos in matter with arbitrary density profiles, 
we adopt the method of series expansion. We suppose that neutrinos are produced as v e and 
then propagate through a given effective potential V(t), this gives the initial values p(0) = 1 
and p(0) = 0. Series expanding the effective potential V(t) and the quantity pit), we obtain 
the following expressions 



p + Gp = 



(10) 



7] 




(11) 



oo 




(12) 



n=0 



oo 




(13) 



n=0 
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where the coefficients V n (n = 0, 1, . . .) define the effective potential and where we wish to 
compute the coefficients p n (n — 0, 1, . . .). By using the relation between / and V, we obtain 

00 A 2 

/(*) = E f» e > ^ = cos ( 2 ^ - v " ( 14 ) 

n=0 

oo n 

F(t) = J2 F nt H , F n = Y,fkfn~k. (15) 
n=0 k=0 

Inserting the above expressions into Eq. (jOJ) and identifying terms of the same order in t 
gives the relation 



F n G = J2( s + 2)(s + l)(n- s + 2)(n- s + l)p s+2 p n - s+2 + 

s=0 
n 

J2 [2G(s + 2)(s + l)p s+2 Pn- s + GW»-J + 

n s 

Fn - S °PkPs-k + {k + l){s-k+ l)p 

k+lPa-k+l- 

(16) 

s=0 k=0 

For n = with the given initial conditions, Eq. (fTBj) is a second order equation in p 2 with 
p 2 = —G/2 as a double root. This corresponds well to the fact that at t — 0, the right-hand 
side of Eq. (JHJ) vanishes for the given initial conditions and we are left with the equation 
p(0) = —Gp(0). For n — 1, Eq. (j!6)l is trivially fulfilled (given the assumed initial conditions, 
terms with p n +2 will appear with the prefactor Gpo + 2p 2 only), while the solution to the 
equation for n = 2 is simply p 3 = 0. 

Also the solution for n = 3 is now trivially fulfilled, since the terms including p n+ i also 
cancel for n > 3. For n — 4, the equation is a second order equation in p 4 with the solutions 

G 2 g(g + Fq) 
P4 = — and p 4 = — . (17) 

Of these two solutions, only the latter will be a solution to our problem, this is easily checked 
by inserting the known solution in the case of constant effective potential from Eq. (fTTj) . 

For n > 5, Eq. (|16J) is now linear in p n . For n > 6, we obtain a solution for p n in terms 
of lower order p k , G, and F s , where k < n and s < n — 4. This expression is the following 
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recurrence relation 



Pn 



G(n 2 -3n + 2)F 

n-2 



n-2 



F 1^2( S + l )( n - S)p s+1 p n - s + 
s=l 

n-3 

G(G + Fq) ^PsPn-s + ^0 + l)(n-S + l)p s +lPn-s+l + 

s=2 s=3 
n—A rt— 1 s 

2G^(S + 2)(s + l)p n _ s p s+2 + G^Fn^s^PkPs-k + 

s=2 s=4 k=0 

n— 2 s 

-k+l + 

s=4 k=0 

n— 3 

y^jn - s + 2){n - s + l)(s + 2)(s + l)f>n- s +2f> s +2 



s=3 



For the first few coefficients we obtain 



Po 
Pi 

P2 
P3 
P4 

P5 

P6 

P7 

P8 



1, 



G 

"2' 



0. 
G 

GF\ 



^(G + F ), 



G(4G F + 8GF 2 + 4F 3 + g - 36F F 2 ) 
2880F 

G(8GF 2 F! + 8F 3 Fi - Ff + 4F FiF 2 - 48F 2 F 3 ) 



5760F 2 



G 



645120F 3 



{48GF 5 + 16F 6 - 63F* + 16F 4 (3G 2 - 34F 2 ) + 



312F F 2 F 2 + 8F (GF 2 - 30F 2 2 - 48FiF 3 ) + 
4F 3 [4(G 3 - 34GF 2 + 240F 4 ) - 53Fi]}. 



(18) 



(19) 



As can be observed by setting F k = 0, the solution for F — 0, i.e., at the MSW resonance, 
is just the series expansion for p = cos(gt), which is clearly as expected. 
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IV. CONVERGENCE OF THE SOLUTION 



In order to test our solution, we perform a number of numerical tests. First of all, we 
give an overview of how we approximate the electron number density (i.e., in principle, the 
effective potential) by a polynomial. Then, we proceed by confirming that our solution really 
converges nicely towards the simple trigonometric function that is the exact solution for a 
constant electron number density. In this case, we also study the convergence of the energy 
dependence of the neutrino oscillation probability P ex (L) for a baseline of L = 3000 km. 
After the constant electron number density case, we investigate the case of a linear effective 
)otential, and finally, we study the case of the Preliminary Reference Earth Model (PREM) 

In the numerical calculations, we have used the mixing angle 9 = 13° and the mass squared 
difference Am 2 = 2 x 10~ 3 eV 2 



19( . The value of 9 approximately corresponds to the upper 
limit on the leptonic mixing angle #13 from the CHOOZ experiment with Am 2 = 2 x 10~ 3 
eV 2 |5] . The reason to use this particular choice of parameters is that #13 and the large mass 
squared difference give the main effects to neutrino oscillations from v e into other flavors for 
the baselines and energies we have studied (for L = 3000 km and E = 1 GeV, the neutrino 
oscillations governed by the small mass squared difference contribute with an approximate 
addition of 0.05 to the neutrino oscillation probability Pp X , for shorter baselines and higher 
energies, this effect decreases), see for example Refs. |20[|. The reason for using the upper 
bound value for the mixing angle 9 and not some smaller value is that we wish to study 
the behavior of our solution rather than to make any precise predictions about the neutrino 
oscillation probabilities. 



A. Series expansion of the effective potential 

In order to use the series solution, which was obtained in the previous section, we will 
need the coefficients V n . In general, for a given baseline length L, the effective potential 
V(t) can be expanded in terms of Legendre polynomials, i.e., 

1 s^Tl O 1 1 /* 1 

V{x) = Y,CnPn{x), P n^) = ^T^[^-m, C n = / V(x)P n (x)dx, (20) 

n=0 ' 



s 



where x = 2t/L — 1. For numerical treatments, we cannot use the entire expansion in Leg- 
endre polynomials because of finite computer memory and finite computer time. However, 
if we assume that the coefficients c n are negligible for n > N, where N is some integer, then 
we have a polynomial approximation 



of the effective potential. Clearly, given any polynomial V(t), it is a trivial matter to extract 
the coefficients V n . This approach turns out to be quite handy in the case of the PREM 
profile, which is discussed below. 

B. Constant matter density 

The first case we study numerically is the case of a constant effective potential. We use 
the baseline length L = 3000 km and the electron number density N e = iV e)Core /3, where 
Kore = V2GFN e;COTe ~ 5.6 x 10~ 19 MeV corresponds to a matter density of about 13 g/cm 3 , 
which is the maximum matter density in the Earth's core Q|. In this case, the coefficients 
V n are easily obtained as V = V(t) and V n = for n > 0. Since the exact solution to 
this problem is known [t|, we focus on the convergence of our solution for P ex , both in the 
energy spectrum and the time evolution. The numerical results are shown in Fig. ^ In this 
figure, we can observe that approximately 20 terms are needed to reconstruct one period 
of oscillation and that the convergence is indeed the same as for a simple trigonometric 
function. 

C. Linearly varying matter density 

Now, we turn our interest towards the case of a linearly varying effective potential. In 
particular, we study a baseline of L = 3000 km, where the electron number density is given 
by N e (t) = N ecorc (t/L + l)/4. As in the case of constant effective potential, it is again 
easy to obtain the coefficients V n from our equation for V(t). Performing the numerical 
calculations result in Fig. |2l In this figure, we have excluded the plot for our solution in 
the shaded region, which roughly corresponds to an energy equal to the resonance energy 
of the effective potential V = Vq, where the solution breaks down numerically. The reason 
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(21) 
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FIG. 1: The neutrino oscillation probability P ex as a function of energy and time, respectively. 
Upper panel: The convergence of the energy spectrum given by our series expansion for a constant 
electron number density profile with N e = iV e)Core /3. Lower panel: The convergence of the series 
expansion for E = 1.2 GeV, corresponding to the bold line in the energy spectrum. The solid 
curves correspond to the exact solutions and the dashed curves correspond to the series expansion. 
The numbers correspond to the number of terms used in the series expansion. The solid vertical 
line in the upper panel corresponds to the energy used for the lower panel. 



for this breakdown can be found in Eq. ()18j) . where we repeatedly divide by Fq. For the 
resonance energy corresponding to V = Vq, we have F ~ 0, which leads to large absolute 
values of numbers that should add up to a number between zero and one. Due to finite 
machine precision, we have numerical errors as a result. 

Apart from neutrino energies near the resonance energy, we can observe that we again 
obtain a nice convergence of both the energy spectrum and the time evolution, where we 
reproduce one full oscillation by approximately 20 terms of our series expansion. It should 
be pointed out that this case of linearly varying matter density has no known application 
to experiments and only serves as an illustrative example. 
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FIG. 2: The energy spectrum of the neutrino oscillation probability P ex for a linear profile (upper 
panel) with N e = N ejCore (l + t/L)/4 along with the convergence of the solution for E = 1.2 GeV 
(lower panel). In the energy spectra, the dotted curve corresponds to the numerical solution for 
the given profile, the dashed curve corresponds to our series solution, where we have included the 
first 35 terms, and the dash-dotted curve corresponds to an approximation of constant matter 
density. The solid vertical line corresponds to the energy used in the lower panel and the series 
solution is not plotted in the shaded region where it breaks down numerically. In the time evolution 
plot, the solid curve corresponds to the numerical solution, the dotted curve corresponds to the 
approximation of constant matter density, the dashed curves correspond to our series solution for 
different numbers of included terms, and the numbers correspond to the number of terms used for 
each of these curves. 

D. PREM profile 

For the PREM electron number density profile, which is the interesting profile in, for 
example, long-baseline neutrino oscillation experiments, we use the expansion in Legendre 
polynomials and truncate the series using N = 2 for definiteness. In effect, this corresponds 
to projecting the function V(t), which is an element in the vector space of real functions on 
the interval [0, L], onto the subspace of second order polynomial functions on [0, L], using 
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FIG. 3: The matter density profiles for different baseline lengths L according to the PREM. The 
solid curves are the exact profiles, the dashed curves are profiles approximated by a second order 
polynomial, and the dotted lines are the average matter densities, i.e., the matter density of 
approximations using constant electron number density. For L = 250 km and 750 km the exact 
profiles and the approximations using a second order polynomial are practically the same and as a 
result they are not distinguishable in the figure. 



the inner product 

(f,g)= [ f(x)g(x)dx. (22) 

JO 

In Fig. El we plot the electron number density profiles for the baseline lengths L = 250 km, 
750 km, 3000 km, and 5000 km along with the second order polynomial approximations and 
the constant approximations. In Fig. 01 we plot the energy spectra and time evolution at 
E = 1.2 GeV for a baseline of L = 3000 km for the PREM profile. Again, our solution is 
not plotted in the shaded region in which it breaks down numerically for the same reasons 
as discussed previously. In this case, it is apparent that if the number of terms used in the 
series expansion is large enough, our solution is a significant improvement from the constant 
matter density approximation. 

Clearly, the approximation of using a second order polynomial for the electron number 
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FIG. 4: The energy spectrum for the neutrino oscillation probability P ex using 35 terms of our series 
solution (upper panel) and the convergence of the time evolution for E = 1.2 GeV (lower panel). 
Here we assume a baseline of L = 3000 km and using the PREM profile for the electron number 
density. In the energy spectrum, the dashed curve corresponds to the numerical solution using 
the PREM profile, the dotted curve corresponds to our series solution, and the dash-dotted curve 
corresponds to the approximation using constant matter density. The solid vertical line corresponds 
to the energy used in the lower panel and the series solution is not plotted in the shaded region 
where it breaks down numerically. In the time evolution plot, the solid curve corresponds to the 
numerical solution, the dotted curve to the solution for the constant matter density approximation, 
the dashed curves correspond to our series solution using different numbers of terms, and the 
numbers correspond to the number of terms used for each of these curves. 



density gives a very good reproduction of the numerical solution (which uses the profiles 
obtained from the PREM). As can be seen in the time evolution plot, the error made is 
barely noticeable until approximately one and a half oscillations, i.e., for lower neutrino 
energies if the baseline length L is kept fixed. This is in good agreement with the results 
obtained in Ref. j^jj, where the effective potential is expanded in a Fourier series, as well 
as Ref. [22], which shows that details of the effective potential that are smaller than the 
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FIG. 5: The energy spectra of the neutrino oscillation probability P ex for different baselines using 
the PREM profile for the electron number density. Again, we have used 35 terms from our series 
solution for each of these energy spectra and the dashed curves correspond to the numerical solu- 
tions, the dotted curves correspond to our series solution, and the dash-dotted curves correspond 
to the constant density approximations. The series solution is not plotted in the shaded regions 
where it breaks down numerically. 



oscillation length cannot be resolved by neutrino oscillations. As noticed in both of the 
earlier cases, about 20 terms are needed in the series expansion in order to reproduce one 
full oscillation. 

For the PREM profile, we are also interested in a number of other baseline lengths. In 
particular, in Fig. |5J we plot the energy spectra for the baseline lengths L = 250 km, 750 
km, 3000 km, and 5000 km. For the baseline lengths L = 250 km and 750 km, there is 
no noticeable difference between the numerical solution, our exact solution, and the ap- 
proximation using constant electron number density. This is to be expected as the electron 
number density does not vary significantly for these baseline lengths (see Fig. EJ). However, 
for both L = 3000 km and 5000 km, we do observe a difference between the constant electron 
number density approximations and the other two solutions. Again, we can conclude that 
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the approximation with a second order polynomial for the electron number density agrees 
remarkably well with the numerical calculation using the profile obtained directly from the 
PREM. Note that the comparison of the energy spectra for different matter density profiles 
and baselines has been studied before 

For baselines longer than L = 5000 km, the region near the resonance for V = Vq tends 
to expand and ruin the numerical convergence of our solution. Also, this region expands 
if we include more terms of the series expansion. For the baseline lengths below L = 5000 
km, this can be somewhat compensated by using fewer terms of the series expansion for 
high energies to avoid the numerical cancellation effects and more terms for lower energies 
to obtain a nice convergence. 



V. OTHER METHODS OF SOLVING THE NEUTRINO FLAVOR EVOLUTION 



In general, there have been numerous methods on how to solve the problem of neutrino 



flavor evolution in matter 



H n, 12, Q, 14, 0, 0, Q. 



obvious formal solution using a time-ordered exponential, i.e., 



W))=T 



exp — i 



H{t)&t 



K0)> 



First of all, there is the 



(23) 



This solution is exact, but it is not very helpful in actual calculations due to the nature of 
the time-ordered exponential. A way of solving this problem is to use a discretization of 
the effective potential The effective potential is then divided into a finite number 

of layers with constant effective potentials {i.e., constant electron number density), which 
approximate a given effective potential. Clearly, when the number of layers goes to infinity, 
one regains the time-ordered exponential. 

There is also the pure numerical approach to the problem, where the neutrino flavor evo- 
lution is easily solved numerically for any effective potential. While this gives the possibility 
of actually calculating numerical values for the neutrino oscillation probabilities, it does not 
offer any insight into how these probabilities vary with different parameters. 



In addition, the neutrino flavor evolution can be analytical 
effective potentials. Examples are the constant (two flavors 
linearly 12] and exponentially varying effective potentials. 



solved for some specific 
and three flavors |ll|). 



Moreover, a widely used solution is the adiabatic solution 



14. 



151 ] . where the effective 
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potential is assumed to change slowly, so that there are no transitions between different 
matter eigenstates of the full Hamiltonian. This approximation can be derived by, for 
example, using the Wentzel-Kramers-Brillouin (WKB) method Q|, where also higher order 
corrections due to non-adiabatic transitions can be calculated. 

There have also been earlier efforts to write the neutrino evolution equations as ordinary 
non-linear differential equations, see for example Ref. However, such equations have 

generally been complex differential equations for the probability amplitudes and not, as in 
the present case, real differential equations for the probabilities. Also, the solutions in such 
cases have been made for special cases of the effective potential and not as series solutions 
valid for all effective potentials. 

Lately, approximate solutions, valid when the effective potential V <C Am 2 /(2E), have 



been presented 



16J] and applied to oscillations of solar neutrinos and the solar neutrino day- 



night effect. However, these solutions are not valid for the baseline lengths and neutrino 
energies which we have treated numerically. 



VI. SUMMARY AND CONCLUSIONS 



We have shown that solving the general problem of two flavor neutrino oscillations with 
an arbitrary effective potential is equivalent with solving the non-linear ordinary differential 
equation 

(p + Gp) 2 = F(t)lG(l-p 2 )-p 2 }, (24) 

where G = \Am 2 /(2E)} 2 sin 2 (2fl), F(t) = [Am 2 cos(2#)/(2£) - V(t)] 2 , and the neutrino 
oscillation probability P ex is given by 

P ex = \{l-P). (25) 

We have presented an exact solution [see Eqs. (fTHj) and (JTHJ)] to this equation by adopting 
the method of series expansion of both the solution and the effective potential V(t) and 
demonstrated the numerical convergence of this solution for a number of different effective 
potentials. In all cases investigated, about 20 terms in the series expansion are required 
to reproduce one period of oscillation. We have also seen that for the neutrino energies 
and baselines considered, the energy spectra of the neutrino oscillation probability is well 
reproduced by approximating the effective potential by a second order polynomial. 
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